home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / lpc / stabilization / xfactor.f < prev    next >
Text File  |  1990-02-03  |  6KB  |  227 lines
  1.       subroutine factor(b,k4,rootr,rooti,kinsid,kprint,eps)
  2. c       sets up problem, calls dproot, 
  3. c       and checks residual values at roots  
  4.       implicit double precision (a-h,o-z)
  5.       double complex z,res,jay
  6.       double precision b(1),rootr(1),rooti(1),coe(37)
  7.       jay=(0.d0,1.d0)
  8.       pi=4.d0*datan2(1.d0,1.d0)
  9.       do 550 i=1,k4
  10. 550   coe(i)=b(i)
  11.       k4m=k4-1
  12.       call dproot(k4m,coe,rootr,rooti,kerr,kprint,eps)
  13. c      write(0,600)kerr
  14. 600   format(' return from dproot with kerr=',i5)
  15.       if(kerr.gt.0)stop
  16.       kinsid=0
  17.       resmax=0.d0
  18.       rmax=0.d0
  19.       rmin=2.d0**(32)
  20.       dist=2.d0**(32)
  21.       amax=2.d0**(32)
  22.       r2=amax**(1.d0/k4)
  23.       do 701 j=1,k4m
  24.       z=rootr(j)+jay*rooti(j)
  25.       r=dsqrt(rootr(j)**2+rooti(j)**2)
  26. c         skip residue calculation if root is too big
  27.       if(r.gt.r2)goto713
  28.       res=b(k4)
  29.       do 705 k=2,k4
  30. 705   res=res*z+b(k4-k+1)
  31.       partr=res
  32.       parti=-jay*res
  33.       resmag=dsqrt(partr**2+parti**2)
  34.       if(resmax.le.resmag)resmax=resmag
  35. 713   if(rmax.lt.r)rmax=r
  36.       if(rmin.gt.r)rmin=r
  37.       if(r.lt.1.d0)kinsid=kinsid+1
  38.       distr=dabs(r-1.d0)
  39.       if(dist.gt.distr)dist=distr
  40. 701   continue
  41. c     write(0,703)resmax
  42. c     write(0,704)rmax,rmin,dist
  43. c703   format('resmax= ',d20.10)
  44. c704   format('rmax= ',d20.10/'rmin= ',d20.10/'dist=',d20.10)
  45.       return
  46.       end
  47.  
  48.       subroutine dproot(mm,a,rootr,rooti,kerr,kprint,eps)
  49. c        mm=degree of polynomial
  50. c        a=coefficient array, lowest to highest degree
  51. c        kprint=1 for full printing
  52. c        kerr=0 is normal return
  53.       implicit double precision (a-h,o-z)
  54.       double complex b(37),c(37),p,pp,z,w
  55.       double complex bb(37),cc(37),jay
  56.       double precision a(1),rootr(1),rooti(1)
  57.       double precision save(37)
  58.       jay=(0.d0,1.d0)
  59.       mmp=mm+1
  60.       m=mm
  61.       mp=mmp
  62.       do 700 i=1,mp
  63. 700   save(i)=a(i)
  64. c       kount is number of iterations so far
  65.       kount=0
  66. c       kmax is maximum total number of iterations allowed
  67.       kmax =20*m
  68. c       newst is number of re-starts
  69.       newst=0
  70. c       ktrym is number of attempted iterations before re-starting
  71.       ktrym=20
  72. c       kpolm is number of attempted iterations before polishing is stopped
  73.       kpolm=20
  74. c       amax is the largest number we allow
  75.       amax=2.d0**(32)
  76.       amin=1.d0/amax
  77. c       rr1 and rr2 are radii within which we work for polishing
  78.       rr1=amin**(1.d0/m)
  79.       rr2=amax**(1.d0/m)
  80. c     eps is the tolerance for convergence
  81.       sqteps=dsqrt(eps)
  82. c        main loop; m is current degree
  83. 10    if(m.le.0)goto200
  84. c        new z, a point on the unit circle
  85.       rkount=kount
  86.       z=dcos(rkount)+jay*dsin(rkount)
  87.       ktry=0
  88. c       r1 and r2 are boundaries of an expanding annulus within which we work
  89.       r1=amin**(1.d0/m)
  90.       r2=amax**(1.d0/m)
  91. c        inside loop
  92. 20    partr=z
  93.       parti=-jay*z
  94.       size=dsqrt(partr**2+parti**2)
  95.       if(size.lt.r1 .or. size.gt.r2)goto300
  96.       if(ktry.ge.ktrym)goto300
  97.       ktry=ktry+1
  98.       if(kount.ge.kmax)goto400
  99.       kount=kount+1
  100. c        get value of polynomial at z, synthetic division
  101.       b(mp)=a(mp)
  102.       do 30 j=1,m
  103.       k=m-j+1
  104. 30    b(k)=z*b(k+1)+a(k)
  105.       p=b(1)
  106.       partr=p
  107.       parti=-jay*p
  108.       if(dsqrt(partr**2+parti**2).gt.amax)goto300
  109. c        get value of derivative at z, synthetic division
  110.       c(mp)=b(mp)
  111.       mdec=m-1
  112.       do 60 j=1,mdec
  113.       k=m-j+1
  114. 60    c(k)=z*c(k+1)+b(k)
  115.       pp=c(2)
  116.       partr=pp
  117.       parti=-jay*pp
  118.       if(dsqrt(partr**2+parti**2).lt.amin)goto300
  119. c        test for convergence
  120.       partr=p
  121.       parti=-jay*p
  122.       size=dsqrt(partr**2+parti**2)
  123.       if(size.gt.eps)goto775
  124.       nroot=mm-m+1
  125.       goto500
  126. 775   continue
  127.       z=z-p/pp
  128.       goto20
  129. c        end of main loop
  130.  
  131. c        normal return
  132. 200   kerr=0
  133.       goto600
  134.  
  135. c        new start
  136. 300   rkount=kount
  137.       z=dcos(rkount)+jay*dsin(rkount)
  138.       ktry=0
  139.       newst=newst+1
  140.       goto20
  141.  
  142. c        too many iterations
  143. 400   kerr=400
  144.       goto600
  145.  
  146. c        root z located
  147. c        polish z to get w
  148. 500   w=z
  149.       kpol=0
  150. 510   partr=w
  151.       parti=-jay*w
  152.       size=dsqrt(partr**2+parti**2)
  153. c       give up polishing if w is outside annulus
  154.       if(size.lt.rr1 .or. size.gt.rr2)goto501
  155. c       give up polishing if kpol>=kpolm
  156.       if(kpol.ge.kpolm)goto501
  157.       kpol=kpol+1
  158.       if(kount.ge.kmax)goto400
  159.       kount=kount+1
  160.       bb(mmp)=save(mmp)
  161.       do 530 j=1,mm
  162.       k=mm-j+1
  163. 530   bb(k)=w*bb(k+1)+save(k)
  164.       p=bb(1)
  165.       partr=p
  166.       parti=-jay*p
  167.       if(dsqrt(partr**2+parti**2).gt.amax)goto300
  168.       cc(mmp)=bb(mmp)
  169.       mdec=mm-1
  170.       do 560 j=1,mdec
  171.       k=mm-j+1
  172. 560   cc(k)=w*cc(k+1)+bb(k)
  173.       pp=cc(2)
  174.       partr=pp
  175.       parti=-jay*pp
  176.       if(dsqrt(partr**2+parti**2).lt.amin)goto300
  177.       partr=p
  178.       parti=-jay*p
  179.       size=dsqrt(partr**2+parti**2)
  180. c       test for convergence of polishing
  181.       if(size.le.eps)goto501
  182.       w=w-p/pp
  183.       goto510
  184.  
  185. c        deflate
  186. 501   b(mp)=a(mp)
  187.       do 830 j=1,m
  188.       k=m-j+1
  189. 830   b(k)=z*b(k+1)+a(k)
  190.       p=b(1)
  191.       rootr(m)=w
  192.       rooti(m)=-jay*w
  193.       m=m-1
  194.       mp=mp-1
  195.       parti=-jay*w
  196.       if(dabs(parti).gt.sqteps)goto140
  197. c        real root
  198.       rooti(m+1)=0.d0
  199.       do 100 j=1,mp
  200. 100   a(j)=b(j+1)
  201.       goto10
  202. c        complex root
  203. 140   partr=z
  204.       parti=-jay*z
  205.       z=partr-jay*parti
  206.       c(mp)=b(mp+1)
  207.       do 110 j=1,m
  208.       k=m-j+1
  209. 110   c(k)=z*c(k+1)+b(k+1)
  210.       rootr(m)=w
  211.       rooti(m)=-(-jay*w)
  212.       m=m-1
  213.       mp=mp-1
  214.       do 130 j=1,mp
  215. 130   a(j)=c(j+1)
  216.       goto10
  217. c        report and return
  218. 600   real1=kount
  219.       real2=mm
  220.       temp=real1/real2
  221. c     write(0,150)kount,temp
  222. 150   format(' kount=',i10,' kount/root=',f15.5)
  223. c     write(0,151)newst,kerr
  224. 151   format(' new starts=',i10,' kerr=',i10)
  225.       return
  226.       end
  227.